Education Statistics in Oregon - Alis volat propriis

Welcome to my notebook covering the use of the tidycensus R package for extracting and evaluating census data in the state of Oregon. In part A We are interested in uncovering the percentages of the population with graduate level degrees at a county level. Part B will continue to focus on education data, Stay tuned to find out more!

Import Libraries

Start by importing the list of libraries that will assist us in connecting with the Census Bureau’s APIs and visualizing the desired data.

library(tidycensus) # connect with the census bureau API
library(tidyverse) # data wrangling
library(mapview) # spatial data visualizer
library(plotly) # graph / chart visualizer
library(ggiraph) # graph / chart visualizer
library(ggplot2) # graph / chart visualizer
library(DT) # fancy tables
library(scales) # scales map data to aesthetics
library(sf) # spatial data package
library(tigris) # erase water
library(cowplot) # plot_grid()

Part A

---

Data Source

We will use the get_acs() function from the tidycensus package to pull data from the census bureau. The function defaults to the 2017-2021 5-year ACS dataset. Which is necessary for us, since there are a list of counties in Oregon with a population of less than 65,000.

The data we would like to visualize is the percent of the population that have a graduate degree at the county level in the state of Oregon.

# use the get_acs() function for importing data
or_grads <- get_acs(
  geography = "county", # county level geography
  variables = "DP02_0066P", # variable for percent of the population that have a graduate degrees
  state = "OR", 
  year = 2021,
  progress_bar = FALSE
)

Data Manipulation & DataTables

An important aspect of the tidycensus package is that is presents the data to us in a tidy format. The data is pre-cleaned and takes little wrangling to answer questions.

# filter and arrange results
or_grads_table <- or_grads %>%
  separate(NAME, into = c("county", "state"), sep = ", ") %>% # create county and state columns
  arrange(-estimate) %>% # sort by estimate descending
  select(county, state, estimate, moe) # select important columns

By using the datatable() function from the DT library we can filter our results to see the desired information.

# function from the DT library to create a datatable
datatable(or_grads_table)

Figure 1 above reveals the data imported from the get_acs() function. Results conclude that the top 3 counties with graduate degrees is Benton, Multnomah, and Washington. Bottom 3 counties are Malheur, Gilliam, and Morrow.

---

Data Visualization & Analysis

Now that our data is imported and in a proper format. We can begin using ggplot to visualize the percentage of the oregon population that graduate have degrees.

# use ggplot with estimate as the x axis and county as the y axis
or_plot <- ggplot(or_grads_table, aes(x = estimate, 
                                      y = reorder(county, estimate))) + # sort by counties high to low
  geom_point(color="navy", size = 3) + # state color navy blue for point data
  scale_x_continuous(labels = function(x) paste0(x, '%')) + # concatenate '%' to estimate data
  scale_y_discrete(labels = function(x) str_remove(x, " County")) + # strip 'county' from the name column
  labs(title = "% Population with graduate degrees, 2017-2021 ACS", # use appropriate title
       subtitle = "Counties in Oregon", # use appropriate subtitle
       caption = "Data acquired with R and tidycensus",
       x = "ACS estimate Percentage",
       y = "") +
  theme_minimal(base_size = 12) # remove tick marks 
  
or_plot # print plot

Figure 2 above plots the counties by estimate value in a descending order. We get a full picture of our results without having to filter through multiple pages in a table.

---

All tables in the tidycensus package come with a margin of error column to give us an understanding of accuracy in the dataset. Typically, census data from geographies of larger population have a lower margin of error.

Lets plot the MOE variable with our estimate variable below.

# create plot revealing moe and estimate variables
or_plot_errorbar <- ggplot(or_grads, aes(x = estimate, y = reorder(NAME, estimate))) + # sort by counties high to low
  geom_errorbar(aes(xmin = estimate - moe, xmax = estimate + moe), 
                width = 0.5, linewidth = 0.5) +
  geom_point(color = "gold", size = 3) + # state color gold for point data
  scale_x_continuous(labels = function(x) paste0(x, '%')) + # concatenate '%' to estimate data
  scale_y_discrete(labels = function(x) str_remove(x, " County, Oregon|, Oregon")) + # strip 'county' from the name column
  labs(title = "% Population with graduate degrees, 2017-2021 ACS", # use appropriate title
       subtitle = "Counties in Oregon", # use appropriate subtitle
       caption = "Data acquired with R and tidycensus. Error bars represent margin of error around estimates.",
       "ACS estimate",
      y = "",
      x = "Percentage") +
  theme_minimal(base_size = 12) # remove tick marks 
  
or_plot_errorbar # print plot

Figure 3 above reveals the margin of error associated with the estimate data.

---

Interative Plots

Using ggplotly() I’m going to recreate the plot above but leveraging the interactivity of the plotly library.

# use ggplotly() and add a tooltip, and extend the margins of the plot
ggplotly(or_plot_errorbar, tooltip = "x")%>% 
  layout(margin = list(l = 50, r = 50, b = 50, t = 50))

Figure 4 Hover over the points in the plot above to see the interactive tootip.

---

Part B

---

Continuing with the theme of education, lets import more data through the tidycensus package.

What percentage of people with bachelor’s degrees rent or own their home in Washington and Multnomah County, Oregon at the tract level?

These two counties together cover the most populus city of Oregon. You guessed it, Portland.

To find variables that will assist us in our analysis, we will use the load_variables() function to filter through the ACS data.

# use the load_variables function() with the year and ACS 5 year parameter
vars <- load_variables(2021, "acs5")
# load vars into a datatable
datatable(vars)

Figure 5 By using the search feature in the datatable, you see examples of available results. For example search “bachelor’s” to find variable codes moving forward.

---

Import Data

We will need to pull in two datasets and eventually merge the results into one table.

“B25013_006” is the population of people with bachelor’s degrees that own their home.

“B25013_011” is the population of people with bachelor’s degrees that rent their home.

# import data at the tract level for home owners
bach_owners <- erase_water(get_acs(
  geography = "tract",
  variables = "B25013_006",
  state = "OR",
  county = c("Washington", "Multnomah"), # select multiple counties
  geometry = TRUE, # get multipolygon data of Oregon tracts
  progress_bar = FALSE
))

# import data at the tract level for renters
bach_renters <- get_acs(
  geography = "tract",
  variables = "B25013_011",
  state = "OR",
  county = c("Washington", "Multnomah"), # select multiple counties
  geometry = FALSE, # geometry is false because we are getting our multipolygon data in the owners dataset. 
  progress_bar = FALSE
  )

Data Wrangling

Use functions of the tidyverse to rename, select, and mutuate columns to fit our analysis. Since the estimate data originally comes as a count, we will sum up renters and owners to a create a total column and a percentage column for each.

# clarify the name of the estimate & moe columns before joining
bach_owners <- bach_owners %>%
  rename(estimate_owners = estimate, moe_owners = moe)

# clarify the name of the estimate & moe columns before joining
bach_renters <- bach_renters %>%
  rename(estimate_renters = estimate, moe_renters = moe) %>%
  select(GEOID, estimate_renters, moe_renters) # drop unnecessary values

# join dataframes on GEOID, and create percent columns of each category
bach_housing <- left_join(bach_owners, bach_renters, by='GEOID') %>%
  select(GEOID, NAME, estimate_owners, estimate_renters, moe_owners, moe_renters, geometry) %>%
  mutate(estimate_total = estimate_renters + estimate_owners) %>% # create a total value
  mutate(percent_renters = (estimate_renters/estimate_total)*100) %>% # create a % of renters
  mutate(percent_owners = (estimate_owners/estimate_total)*100) %>% # create a % of owners
  mutate(across(where(is.numeric), ~ round(.x, digits = 2)))

Error in exists(cacheKey, where = .rs.WorkingDataEnv, inherits = FALSE) : 
  invalid first argument
Error in assign(cacheKey, frame, .rs.CachedDataEnv) : 
  attempt to use zero-length variable name

Data Visualization & Analysis

Now that the data is wrangled we can view it in a table before mapping and plotting the desired variables.

# select percent renters and owners
bach_housing_table <- bach_housing %>%
  select(NAME, percent_renters, percent_owners) %>%
  arrange(-percent_renters) # arrange in descending order

# preview in datatable
datatable(bach_housing_table)

Figure 6 As suspected in Portland there are a high number of renters with bachelor’s degrees. Multiple tracts reveal areas with 100% of our particular study rent their homes.

We are going to start with the mapview package to get a glimplse into the data.

# use the percent_renters variable to display in the mapview
mapview(bach_housing,
        layer.name = "% Renters<br/>2017-2021 ACS", # change layer name
        zcol = "percent_renters") # select a variable to display

Figure 7 Mapview is programmtically a simple way to display insightful and interactive data. We can now easiy visualize that the city center has the highest number of renters versus owners and shifts as you move into the more rural areas.

Now that we have sampled the data in mapview, lets use ggplot to create a set of chloropleth maps for both categories (renting and owning).

# create ggplot for renters
renters <- ggplot(bach_housing, aes(fill = percent_renters)) + # fill with percent renters
  geom_sf() + 
  theme_void() + 
  scale_fill_viridis_c(option = "rocket") + # use rocket color scheme
  labs(title = "Percent of Population w/ Bachelor's Degree Renting by Census tract", # create clear titles and captions
       subtitle = "Washington & Multnomah County, Oregon",
       fill = "% Renting",
       caption = "2017-2021 ACS | tidycensus R package")

# create ggplot for owners
owners <- ggplot(bach_housing, aes(fill = percent_owners)) + 
  geom_sf() + 
  theme_void() + 
  scale_fill_viridis_c(option = "rocket") +  # use rocket color scheme
  labs(title = "Percent of Population w/ Bachelor's Degree That Own Their Home by Census tract",  # create clear titles and captions
       subtitle = "Washington & Multnomah County, Oregon",
       fill = "% Owners",
       caption = "2017-2021 ACS | tidycensus R package")

# use plot_grid from the cowplot library to stack the maps
plot_grid(renters, owners, ncol=1, align = "hv")

Figure 8 The last figure paints a clear picture of the city center vs rural areas in terms of the ability to own a home in the Portland Metro area.

---

---